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We discuss in detail properties of trapped atomic condensates with anisotropic dipole interactions. 
A practical procedure for constructing anisotropic low energy pseudo potentials is proposed and 
justified by the agreement with results of numerical multi-channel calculations. The time dependent 
variational method is adapted to reveal several interesting features observed in numerical solutions 
of condensate wave function. Collective low energy shape oscillations and their stability inside 
electric fields are investigated. Our results shed new light into macroscopic coherence properties of 
interacting quantum degenerate atomic gases. 
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I. INTRODUCTION 



Q\ . The recent success in atomic Bose-Einstein condensation (BEC) has stimulated great research activities 

t-H ' into trapped quantum gases B. To a remarkable degree, a single condensate wave function of the macroscopically 
occupied ground state, described by the nonlinear Schrodinger equation (NLSE) [H, captures all essential features 
of its coherence properties ||. In fact, one of the key diagnosis features for BEC, the reversed aspect ratio of a 
free expanding condensate, is described purely by the condensate wave function In the standard treatment for 

' the condensate wave function of interacting atoms, realistic inter-atomic potential V{R) is often not directly used. 
Instead, a contact pseudo potential, UqS(R), obtained under the so-called shape independent approximation (SIA) |^| 
is used. Such an idealization results in tremendous simplification, yet to date, SIA has worked remarkably well as 
verified by both theoretical calculations and experimental observations |i|, |lo| , |l~i]| . 

Currently available degenerate quantum gases are cold and dilute, the interaction is therefore dominated by s-wave 
collisions, described by a single atomic parameter: a sc , the s-wave scattering length, if the inter-atomic potential is 
■ isotropic and short ranged (decaying fast than —l/R 3 asymptotically). The complete scattering amplitude is then 
isotropic and energy- independent, given by f(k, k') = — 47ra sc for collisions of incoming momentum k state scattering 
into k'. 

One of the attractive features of atomic degenerate gases lies at effective means for control of the atom-atom 
interaction |p^-p^|. Indeed, very recently several groups have successfully implemented Feschbach resonance |l5| , |l6| , 
£j ', thus enabling a control knob on a sc through the changing of an external magnetic field. Other physical mechanisms 
r* ■ also exist for modifying atom-atom interactions, e.g. the shape resonance as proposed in Jl^| . In an external electric 
. £h \ field, inter atomic potential is modified by the addition of an anisotropic (induced) dipole interaction. 

Although anisotropically interacting fermi system has been an important area of study, e.g. liquid 3 He |1| and 
d-wave paired high T c superconductors |fl9|| . Its bosonic counterpart has not been studied in great detail. In particular, 
we are not aware of any systematic approach for constructing an anisotropic pseudo potential 

For bosonic systems, another related topic is the condensate stability. Under the SIA, the scattering length takes a 
positive or negative value, corresponding to repulsive or attractive interactions. When a sc < occurs, self-interaction 
leads to a collapse of BEC in dimensions higher than 1 pp| , thus the resulting condensate is limited by a critical 
number of particles plj . Anisotropic dipole interactions, on the other hand, are more complicated as both attractive 
and repulsive interactions arise along different directions. We note that several recent investigations have studied 
efforts of non-local interactions on condensate stability p^ , p3[ . 

In this paper, we study properties of trapped BEC of atoms with dipole interactions j24j], arising from either 
external electric field (induced) or permanent magnetic moments p5| . We propose a practical method for constructing 
anisotropic pseudo potentials that can also be extended for investigation of polar molecular BEC |2l],^8| . This paper 
is organized as following. We first briefly review the usual pseudo-potential approximation under the SIA. In Sec. 
II we describe and justify in detail a procedure for constructing effective low energy pseudo-potentials of anisotropic 
interactions. In Sec. Ill we provide our numerical procedure for solving the NLSE with anisotropic dipole interactions. 
Particular emphasis is put on the careful treatment of the singular origin of dipole interactions. We also present and 
discuss results from selected numerical calculations. To explain the stability region as well as the interesting aspect 
ratios observed from our numerical calculations, we perform in Sec. IV an analytic time dependent variational 
calculation. We compare the results obtained with direct numerical solutions of NLSE. Finally we conclude. 
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II. FORMULATION 



For N trapped spinless bosonic atoms in a potential Vt(r), the second quantized Hamiltonian is given by 



n = / df¥( 



■^ 2 + V t (r)-» 



+ \ J dft j ^^(^(r'Wir- f')*(r f )*(r), (1) 

where &(r) and w^(r) are atomic (bosonic) annihilation and creation fields. The chemical potential \i guarantees the 
atomic number N = J df^f' (r)W(f) conservation. 

The bare potential V(R) in (Q) needs to be renormalized for a meaningful perturbation calculation The usual 
treatment is based an effective interaction obtained by a resummation of certain classes of interaction diagrams [p9|-|3l| . 
Physically the SIA can be viewed as a valid low energy and low density renormalization scheme, one simply replaces 
the bare potential V(R) by a pseudo potential uqS(R) whose first order Born scattering amplitude reproduces the 
complete scattering amplitude (— a sc ). This gives u = 4Trh 2 a sc /M. 

When an electric field is introduced along the positive z axis, an additional term 

V E (R)=u 2 ^, (2) 

appears in the atom-atom interaction |l7|, where 112 = — 4y/ (7r/5) a(0)a* (0)£ 2 . a(0) is the atomic polarizability, and 
£ denotes the electric field strength. As was shown before [ fi7| , this modification results in a completely new form for 
the low-energy scattering amplitude 

f(k,k') = 4tt J2 t l t\m* m (k)Y, ml (k'), (3) 

lm,l'm' 

with t\™ (£) the reduced T-matrix elements. They are all energy independent and act as generalized scattering 
lengths. The anisotropic Ve causes the dependence on both incident and scattered directions: k and k 1 = R. 
We therefore propose a general (energy-independent) anisotropic pseudo potential constructed according to 

V eS (R) = uo5(R)+ J2 (4) 

li >0,mi 

whose first Born amplitude is 

/Bo r n(M') = -(4^) 2 a sc y * (fc)y 00 (fc') - ^2 E 7; I - 1 (4 7 r) 2 EE^ m '^' m i) y 4(^' m '(fc / ), (5) 

linii Im I'm 1 

with T£»'(h, mi ) = (i) l + l 'TZ\'l l l '^'(l 1 ,m 1 ). Both 
l\'™'(hm x )= {Y Vm ,\Y limi \Y lm ) 



, ,,,„ / (2/ + l)(2^ + l)(2/ 1 + l) ( I V h\ I V in 
[ L) V 4tt I -m m! m x )\ ' ' ' ' 



and 



r°° 1 
l = J dR ^ kR )3^ k ' R ) 



7T , T( lJ f) _ f-l-l' + l l + l' , 3 2 \ 

= ^ + |) ^ ' ~ ? + 2^ J ■ (?) 

can be computed analytically. The 1/R? form in Eq. (|4|) assures all TZ[ to be k = k' independent (easily seen by a 
change of variable to x = kR in the integral). Putting 
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/Born (MO = /(MO, (8) 

i.e. requiring the Bohn amplitude from the pseudo potential Eq. (^) to be the same as the numerically computed 
value f(k, k'), one can solve for all 7/ lTOl (£) from known t\™ (£) |l7],|3^]. This reduces to a set of (under determined) 
linear equations 

"^E l hm A^)Ti™' {h,m x ) = t l f, (9) 

Iiirti 

for all (Im) and (I'm!) with /, /' 7^ 0, and separately a sc (£) — — ). 

Considerable simplification arises further for bosons (fermions) as only even (odd) (I, I') terms are needed to match 
in (^). Figure |l| displays the result of field dependent a sc (£) 41 K atoms in the triplet electron spin state. Note the 
spikes of shape resonances. The Born amplitude for the dipole term Ve is 

/ Born (fc,fc')=«2-^(47r) 2 T 2 £ T^Y^Cm^Ck/), (10) 

lm,l m 

with T 2 ° = -0.023508. T^' = T^ n '(2, 0)/T 2 ° are independent of electric field 5 within the perturbative Born 
approximation as tabulated below. 



TABLE I. = t\'™'(£)/C(£) for small (1,1'). 



(lm),(l'm') 


(00) 


(20) 


(40) 


(60) 


(80) 


(00) 





1 











(20) 


1 


-0.63889 


0.14287 








(40) 





0.14287 


-0.17420 


0.05637 





(60) 








0.05637 


-0.08131 


0.03008 


(80) 











0.03008 


-0.04707 



We find that away from regions of shape resonances to be discussed elsewhere [ j32| , all numerically computed t\™ (£) 
values, large enough to justify their inclusions, are actually all proportional to £ 2 . Thus, we could rewrite f(k,k') as 

f(k,k') = (iir)C(£) £ t l ^'YC m (k)Y Vm ,(k% (11) 

with scaled quantities t\™ = t\™ (£)/*oo(£) now au being constants. We have since computed (numerically) for 
several alkali metal isotopes, our results for t\™ are tabulated below. 



TABLE II. 4m ' for 7 Li. 



(lm),(l'm') 


(00) 


(20) 


(40) 


(60) 


(80) 


(00) 





1.0 


0.0 


0.0 


0.0 


(20) 


1.0 


-0.63 


0.14 


0.0 


0.0 


(40) 


0.0 


0.14 


-0.17 


0.057 


0.0 


(60) 


0.0 


0.0 


0.057 


-0.080 


0.031 


(80) 


0.0 


0.0 


0.0 


0.031 


-0.044 
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TABLE III. 


tZ for 39 K. 






(lm), (l m ) 


fr\r\\ 

(00) 


(20) 


(40) 


/ r>c\\ 

(60) 


(80) 


(00) 





1.0 


0.0 


0.0 


0.0 


(20) 


1.0 


-0.64 


0.14 


0.0 


0.0 


(40) 


0.0 


0.15 


-0.17 


0.085 


0.0 


(60) 


0.0 


0.0 


0.085 


-0.127 


0.047 


(80) 


0.0 


0.0 


0.0 


0.047 


-0.074 







TABLE IV. 4T 


for 41 K. 






(lm),(l'm') 


(00) 


(20) 


(40) 


(60) 


(80) 


(00) 





1.0 


0.0 


0.0 


0.0 


(20) 


1.0 


-0.64 


0.14 


0.0 


0.0 


(40) 


0.0 


0.14 


-0.17 


0.057 


0.0 


(60) 


0.0 


0.0 


0.057 


-0.081 


0.030 


(80) 


0.0 


0.0 


0.0 


0.030 


-0.047 







TABLE V. 47 


for 85 Rb. 






(lm),(l'm') 


(00) 


(20) 


(40) 


(60) 


(80) 


(00) 





1.0 


0.0 


0.0 


0.0 


(20) 


1.0 


-0.64 


0.14 


0.0 


0.0 


(40) 


0.0 


0.14 


-0.17 


0.056 


0.0 


(60) 


0.0 


0.0 


0.056 


-0.081 


0.030 


(80) 


0.0 


0.0 


0.0 


0.030 


-0.047 







TABLE VI. t\™ 


for 87 Rb. 






(lm),(l'm') 


(00) 


(20) 


(40) 


(60) 


(80) 


(00) 





1.0 


0.0 


0.0 


0.0 


(20) 


1.0 


-0.64 


0.15 


0.0 


0.0 


(40) 


0.0 


0.15 


-0.17 


0.056 


0.0 


(60) 


0.0 


0.0 


0.056 


-0.079 


0.032 


(80) 


0.0 


0.0 


0.0 


0.032 


-0.045 



The agreement between the first order Born approximation and the multi-channel scattering calculations is remark- 
able. We estimate the numerical scattering results to be accurate to a few percent (except for 39 K), independent of 
atoms being bosons (even I, V) or fermions (odd I, /'). Only bosonic results are being considered in this paper. This is 
displayed by noticing the agreement between Table [j| (~ a few per cent) with Table |J-|V]]. This interesting observation 
applies for all bosonic alkali triplet states: 7 Li, 39 ' 41 K, and 85 ' 87 Rb, for up to a field strength of 3 x 10 6 (V/cm) 0§2| 
computed by us. Physically, this implies the effect of Ve is perturbative when £ remai ns s mall (in a.u.). For the 
convenience of further discussions, we tabulate polarizabilities of selected atoms in Table VII. 



TABLE VII. Atomic polarizabilities in (a.u.). 





H 


Li 


Na 


K 


Rb 


a(0) 


4.5 


159.2 


162 


292.8 


319.2 
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What is remarkable is the fact that Tq®{£) and t 2 S,(£) also agree in absolute values |[7| except for a slight difference 
(1-6%). They are calculate below and tabulated in VTJT 



20_ 
00 " 



-16tt 
1718 x Ma?Ta 



(12) 



where all overlined quantities are in atomic units (a.u. 



TABLE VIII. Comparison of numerical values of u 2 ^ i {4'k) 2 T^ with — (47r)t§§(£), the cause of the slight difference is 
unclear but within numerical errors. 





Born 


Multichannel 


atom 


u 2l ^(4nfT f(a ) 


-(4tt)^8(£)(oo) 


7 Li 


3.040 x 10 8 ? 2 


3.238 x 10 8 ? 2 


39 R 


5.713 x 10 9 £ 2 


5.699 x 10 9 ? 2 


41 K 


6.006 x 10 9 ? 2 


5.99 x 10 9 £ 2 


85 Rb 


1.486 x 10 10 £ 2 


1.474 x 10 10 £ 2 


87 Rb 


1.495 x 10 10 £ 2 


1.512 x 10 10 £ 2 



An important parameter in our discussion is the ratio between u 2 and Uq. Since the results from first order Born 
approximation and the multichannel calculations are about the same, one can write this ratio as 

c(£) = -^ =1 £ 2 , (13) 

UQ 

where 7 ~ 1.748 x 10~ 17 a 2 M/a sc and £ is in unit of V/cm. The 7 values for selected atoma are tabulated in Table 
IX. 



TABLE IX. Values of parameter 7 for selected atoms with their field free (E = 0) scattering lengths. 





H 


7 Li 


23 Na 


41 K 


87 Rb 


7(xl0" 13 ) 


2.9 x 10~ 3 


-1.1 


2.0 


2.2 


15.0 
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Most atoms possess permanent magnetic dipoles. With alkali metals, the magnetic dipole mainly originates from 
valance electron spin, typically measured in units of Bohr magneton. It is therefore interesting to compare electric 
dipole interactions with magnetic dipole interactions. The dipole interaction strength between atoms of a permanent 
magnetic dipole \l is 



2—2 2 

= 1.331 x l(T 5 7Z 2 e 2 a§, 

with [Ib the unit of Bohr magneton. For induced electric dipoles, the interaction strength is 

2 C 2 — 2"c 2 2 2 

at — ate a n . 



(14) 



(15) 



A typical heavy alkali atom has a ~ 200, thus for which a 1 (hb) magnetic moment corresponds to an effective electric 
field of 3.3 x 10 -5 (a.u.), or 1.7 x 10 5 (V/cm). Atoms with large magnetic moments effectively simulates induced 
dipole interactions at a high value of equivalent electric field Q . 

As can be concluded from comparing listed data in all tables, we can approximate Eq. ([!]) by keeping only the 
l\ = 2, mi = term to achieve a satisfactory level of accuracy. Thus away from shape resonances we take 



V eS (R) = u S(R) + u 2 Y 20 (R)/R 3 , 
where uq — 4 ^ a sc (£) and U2 = —c(£)uq. The Hamiltonian (|l]) then becomes 

h 2 



(16) 



H 



2M 



V 2 + V t (r)-[i 



u 



(17) 



with R = f — ¥ . The Heisenberg equation for ^>(f,t) becomes nonlocal. At zero temperature the condensate wave 
function ip(r,t) = obeys the NLSE 



-|^V 2 + V t (r) - n + u \^(r 7 t)\ 2 + u 2 J dr* 



Y 20 (R) 
i? 3 



(18) 



with tp(f,t) normalized to N (the number of the atom in the condensate). 



III. NUMERICAL STUDIES 



In this section we discuss the ground state properties of trapped condensates based on numerical solutions of NLSE 
(|18|). We start with a detailed analysis of our numerical procedure for handling the non-local dipole interaction 



A. The numerical procedure 



We use steepest descent through a propagation of Eq. ( jig ) in imaginary time (it) to find its ground state wave 
function. With an axial symmetric harmonic trap 



V t (r) = -Mv 2 (\lx 2 + \ 2 yy 2 + A 2 z 2 ), 



(19) 



we rescale Eq. (|iq ) by introducing dimensionless units for length (at = y/h/Mv), energy (hi/), time (2i/v), and wave 
function (yiV/af). We than obtain 



(20) 



G 



with 



ft = -V 2 + (A 2 .x 2 + X 2 y y 2 + A 2 z 2 ) - 2 M + 2{2nf' 2 P 



\mt)\ 2 -c(£) 



(21) 



where P = y/2/nNa sc /a t and ip(r,t) is normalized to 1. 

The ground state is found be starting with an arbitrary random wave funtion, and propagating Eq. ( |2p| ) in t until 
it's stable (apart from its decreasing norm). In practice we chose an appropriate time step At and iterates the Eq. 
( pp| ) according to 



ip(r,t + At) = ip(r,t) - (At)Hip(r,t). 



(22) 



We renormalize ip to 1 after each iteration and adjust At to control the rate of convergence. 

For a cylindrical symmetric trap (A^ = A y = 1, X z = A), the ground state wave function also possesses the cylindrical 
symmetry. Therefore the non-local term simplifies to 



df*Mp',z')f 



. Y W {R) 
R 3 



dz> / dp'K{p',p,z' -z)\^{p',z'tf 



(23) 



with a kernel 



fc(p',P,z' - z) = ~\ - 



n [(P - P') 2 + ~ z) 2 ] 2 [{P + P') 2 + 0' ~ z) 2 p 
[(p 2 -p' 2 ) 2 -2(p 2 +p /2 )(z'-z) 2 -3(z'-zf]E 

w 



(p + p') 2 + {z'-zf 



+ [{p~p') 2 + {z' ~z) 2 ]{z' -z) 2 K 



[{P + P') 2 + {z 1 - zf 



(24) 



where E[.] and K[.] are standard Elliptical integrals. 

We discretize the {p, z) plane into a two-dimensional grid of points such that wave function values at each point 
becomes a matrix. At each time step the matrix elements are updated according to (|22|). The derivatives in the 
Hamiltonian are evaluated by means of finite-difference methods. Typically, the ground state can be sufficiently well 
described using a grid of 100 x 200 points in the range < p < 5 and —5 < z < 5. 

At first sight, one may naively underestimate the complication of Eq. ( po|) due to the non-local interaction term 
in (|2l]). Several other groups have addressed non-local interactions previously pT| , p3[ . There is however, a significant 
numerical challenge with the dipole interaction, which is singular at the origin. To accurately represent its detailed 
structure, an enormously large grid is needed. Although a Fourier transform into momentum space could simplify the 
convolution operation of the nonlocal term. We found it hard to completely avoid the effort of the singularity this 
way by going to a momentum representation with a limited coarse grid |25|] . Physically, this singularity implies the 
presence of two different length scales for Eq. ( pp[ ) . We thus developed a numerical procedure with two overlaying 
grids: a coarse grid for the relatively smooth wave function and a much finer grid for computing the non-local dipole 
interaction kernel. 

The kernel tC(p', p, z' — z) is divergent at r — f, we thus define a cut-off radius R c such that K(p', p,z' — z) = 
whenever \r — r"| < R c . This cut off is chosen to be small enough that negligible errors result from the numerically 
represented kernel. Typically R c ~ 50(ao) taken, which is much smaller than the wave function grid size. The rapid 
varying kernel fC(p', p, z' — z) is treated with a finer grid. Instead of directly integrating over K.(p', p, z' — z)\ip(p', z')\ 2 
on the wave function grid, we first integrate the kernel separately over the fine grid around each of the wave function 
grid point. Such an integration is numerically intensive, but only needs to be performed once for each of the wave 
function grid point as the kernel is determined by the geometry of system. The integrated kernel values on the 
wave function grid remain the same for different traps and different number of atoms. Finally the non-local term is 
approximated by integrating over the wave function grid using the product of integrated kernel values and the wave 
function. 

For a homogeneous distribution of aligned dipoles, the mean dipole interaction vanishes as Y2o(R) averages to zero 
upon integration over df or df'. This property is maintained for our kernel Eq. (p4|) even though we have integrated 
over (</> — </>') first. We have verified this by noting that the integration of K.(p' , p, z' — z) over (p, z) and (p\ z') does 
vanish. 

For cylindrically symmetric traps, the wave function grid as well as the integrated kernel region is as illustrated in 
Fig. [|. An accurate representation of the integrated kernel requires a quadrature operation over a much finer grid for 
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each of the shaded regions surrounding the wave function grid. As is shown in Fig. g, there are three different types of 
shaded regions, labeled as 1, 2, and 3. Both types of 1 and 2 are boundary terms, which are not needed since they are 
respectively at z = ±L Z and p = 0, L p where either the wave function ip(p' > z') varnishes or the integration measure 
J p 'dp' vanishes. Therefore, we need only to compute the integration of kernel over type 3 element by defining a 
much finer grid and use standard numerical quadrature techniques. When (p, z) ^= (p' , z'), the integration reduces to 
a two-dimensional one which can be easily performed. On the other hand, a three-dimensional integration is needed 
when (p, z) = (p', z'). In this case we have to carefully implement the cut-off radius R c . 

Figure |3J compares the kernel with the coarse grained integrated kernel. Note the significantly different vertical 
scale. 



B. Vortex states 

Simple vortex states ]33]-|3l| with quantized circulations can also be considered by writing the condensate wave 
function in the form 

<P(r) = \^{r)\e- m \ (25) 



with (f> the azimuthal angle with respect to z-axis. The corresponding Eq. ( |2C| ) for \ip(r)\ is then modified by the 
addition of n 2 / p 2 to H of Eq. @. 



C. Numerical Results 



1. Ground state wave function 

The ground state properties of trapped condensates with dipole interactions were first discussed by us in pj| ] . Our 
basic findings are: 1) condensates become elongated along the direction of external field while shrank in the orthogonal 
radial direction; 2) for given values of P and A, there exists a maximum cm(P, A) of a threshold field strength beyond 
which condensate collapse occurs. Figure [| displays numerically computed cm for A = 1 at several different P values. 
For comparison we also show the results from variational calculations to be discussed later. 

The condensate collapse is mainly due to attractive dipole interactions along the direction of external field. To 
minimize its total energy, condensed atoms prefer to align along the attractive direction (z-axis), while narrowing its 
width along the radially repulsive direction. The collapsing occurs when the radial width eventually approaches zero 
with increasing external field strength |36| . 

Figure]^ shows a vortex state wave function for n = 1. The effects of dipole interactions are similar to the ground 
state. As c(£) increases, the vortex state will also collapse. Because of the zero density inside the vortex core, cm in 
this case is much larger than for the ground state. 



2. Comparing numerical solutions with TFA 

A useful approximation for the ground state solution of NLSE ([II]) is the Thomas-Fermi approximation (TFA). 
When the interaction between atoms are repulsive, i.e. with a positive scattering length, the condensate is expected to 
increase its size as compared to the single atom ground state in the trap. With more atoms, the larger the condensate 
size, eventually the spatial directives, consequently the kinetic energy term becomes negligible. In this limit, TFA is 
used to find the ground state wave function by neglecting the kinetic energy term in Eq. (^l|) . With a SIA interaction 
term, the solution simply takes the shape of the inverted trap potential Vt( f). The nonlocal dipole interaction term, 
however, prevents a simple analytic solution even with the TFA. 

Typical solutions to Eq. jTsl ) with and without TFA for P = 5000, A = 1, and c(£) — 0.2, 0.6, 0.7 are compared in 
Fig. | 



IV. TIME-DEPENDENT VARIATIONAL ANALYSIS 



The time-dependent variational approach can also be used to analyze solutions of (|18|) |p7 38 1. We start by identi- 
fying a Lagrangian density C 
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C = -h 

2 



u ^ dip*(f) dip(r) 



Of 



+ ^ M \^(r)\ 2 + Vt(r)\m\' 



u 



\m\ 4 



"2 



2 ,rv " ' 2 IT v " 7 i? 3 
The NLSE can then be found from a minimization of the action p8[ 

5 = / £dfdt. 



(26) 



(27) 



To simplify the variational calculation, we restrict ip to a convenient family of trial functions and study the time 
evolution of the parameters that define the family. A natural choice is a Gaussian-like function first used in |B8[ 



fl>(x,y,z,t) = A(t) J] 



: ,-[ri-ri o (t)] 2 /2w*+iria v (t)+ir i 2 T ,(t) 



(28) 



r]—x : y^z 



where A (complex amplitude), w v (width), a v (slope), (3 V (curvature radius) -1 / 2 , and rjo (center of cloud) are vari- 
ational parameters. This approach, pioneered by Perez-Garcia et al. pq ], has since been successfully used for many 
studies of trapped condensates, a more recent application attempted to explain the anomalous behavior in the finite 
temperature excitation experiment [|39|]4C| 1 . 

Our goal here is to find equations governing the variational parameters. To this aim, we insert ( p8| ) into ( p6| ) and 
calculate an effective Lagrangian L by integrating the Lagrangian density over all space coordinates 



L=(£)= / 

J — i 



Cdr. 



(29) 



to arrive at 



L= 



ihN / A* 



N 



E 



A* A 



2h 2 (3 2 v 1 



N 2 



2\f8TT 3 / 2 W x W y W z 



v[ + -M v z \ 2 j (w 2 + 2rj 2 ) + ^ + 

l n 'I 



M 

3cos 2 6» - 1 



277 



2Mw 2 



H 2 a 2 
M 



where we have used atom number conservation 

N = TT 3/2 \A(t)\ 2 w x (t)w y (t)w z (t) = const. 



(30) 



(31) 



At this point, the variation calculation effectively has been reduced to a finite dimensional problem, i.e., to solve 
the Lagrange equations 



dL 



dL 



0, 



dt \d(jj J dqj 
with the notation 

qj = {w X: w y: w z ,A 1 A*,xo,yo,zo,a X: a y ,a z ,p x: Py : p z }. 
We find equations for the center of the condensate 

Vo + -\V 2 »7o = 0, 



(32) 

(33) 
(34) 



and the condensate widths satisfy 
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w„ + A 2 ^ 2 w„ = 



N d 



71 ^ " " " M*w* AV2^/ni dw v 
The rest of the parameters can be obtained from 

and 



_ MW V ^ 

11 2hw v ' 



lfc/^ eXP {-^^} 



3cos 2 6> - 1 



'/<> 



770^ 



It is convenient to introduce new dimensionless variables t — vt and w n = a t u r; . We then arrive at 



dr 2 



u», + \v n = — - P 



a 



1 



V X VyV Z 



1 - 



5^)/«^.{-x:^} 



3cos 2 # - 1 



(35) 



(36) 



(37) 



(38) 



where P = ^2/itNa sc /at- This equation describes the motion of a particle with coordinates (i? x , v y ,v z ) in an effective 
potential 



U(v„ v„v,)= i (A J ,o» + A^-, 2 + Ajof) + I (± + ± + ±j 

L c(£) | ir - exp {_£|l} 



_5_ 

16? 



3cos 2 6> - 1 



(39) 



V X VyV Z 

For a cylindrically symmetric trap with \ x = X y = 1, X z = A, all integrals can be performed analytically to yield 

1 P 

(40) 



d 2 IP 
— « + w = - + — I 1 - c (^)/(«)] 



d 2 IP 

— T n 2 + \ 2 v z = -3 + [1 - c(£)#(k)] , 



(41) 



with Wj; = v y = v, k = v/v z , and 

m 



5tt 



6(1 -k 2 ) 2 



[-4k 4 - 7k 2 + 2 + 9k 4 #(k)] , 
[-2k 4 + 10k 2 + 1 - 9k 2 P(k)] . 



3(1 -k 2 ) 2 

H(n) = tanh _1 Vl — k 2 /%/1 — k 2 . The equilibrium widths are then determined by 



(42) 



1 P 

u o= — + — I 1 - c(£ )/(«o)] , 

1 P 
A 2 ^o= -3- + -5-3- [1 - c(£ )<?(k )] 

y z0 u u z0 



(43) 



Using w 2 o and ko, Eq. (43) can be rewritten as 

KqVzO 



1 



3„3 1 „3„4 



K w z0 ^O^zO 



[1 - c{£ )/(«o)] , 



A i> z o= 



1 



P 



^zO ^O^zO 

In following discussion, we consider only a sc (£) > case, which implies both P > and c(£) > 0. 



[1 - c(£)g(K )] 



(44) 
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A. Equilibrium widths 



First, for simplicity, we assume that our system satisfy Thomas- Fermi limit (?> 1), then we can safely ignore the 



kinetic term and rewrite Eq. (44) for ko in the following form 

4 [1 - c(Off(«o)] = A 2 [1 - c(£)/(«o)] ■ (45) 

This equation can be solved graphically. From Fig. ^, we first note both f{n) and g{n) are monotonically decreasing 
functions bounded between t/5tt/3 and — 2\/57r/3, also the inequality /(«) > g(n) holds for all n > 0. Then for all 
k > A, we have k 2 [1 — c{£)g{n)} > A 2 [1 — c{£ )/(«)], therefore, if k is a solution of Eq. (|45|), it must satisfy kq < A. 
Meanwhile, when c(£ ) = 0, the solution for « in Eq. ( flBf ) is ko = A. This then proves that no matter what the initial 
field- free condensate aspect ratio is, the condensate always become more prolate along the electric field direction, i.e. 
approximately as illustrated in Fig. |^, it expands along the field direction but shrinks in the orthogonal direction. As 
will be discussed in more detail later, the total condensate volume actually shrinks with increasing fields because of 
the attractive dipole interaction. 

This result, first explained by us p4| ] in terms of the minimization of total energy, is different from the conclusion 
reached in Ref. based on a force argument. This interesting feature has also been independently verified by 
numerical solutions based on the FFT algorithm adopted by Goral et. al p5fl . 

We can also rewrite Eq. (Ea) as 



,2 \2 



A = c{£)h{n), (46) 

with h(n) = K 2 g(n) — A 2 /(k). Figure shows its graphical solutions (A = 1) at several different c(£) values. 

First, we note that h(n = 0) = — v5n\ 2 /3. Thus, as long as c{£) < 3/V5n, there will be one and only one root. 
This result may not look absolutely clear from the figure because of plotting constraints, but it can be seen clearly 
form Eq. (f45|). We also find that as c(£ ) increases, there may exist one, two, or zero roots for k. Once ko is known, 
one can easily find solutions vq and v z q, whose stability can also be checked straightforwardly. It turns out that, in 
all our calculations, whenever only one root for k occurs its corresponding solution for vq and v z q is always stable. 
If there are two roots for k, the solution for vq and v z q corresponding to the smaller no root is a saddle point, thus 
always unstable, while the other is always stable. 



We now consider the A dependence of the condensate property. Figure 10 shows the function h(n) at several different 
A values. We see that as A increases, the maximum value of function h[K), h max also increases, and ft. max = at 
A c = 5.170169. For those A values corresponding to a negative h max , no root of k is found if c(£) becomes sufficiently 
large. But when A > A c , at least one root for k exists no matter how large a c(£) is. This means the condensate 
cannot simply collapse even at these very large electric field strengths. Physically this implies the increased stability 
of condensate inside an electric field with increasing values of A. A condensate of a pancake shape is more stable than 
one with a cigar shape. This can be simply understood from the following argument; The collapse of a condensate 
under electric fields is due to the alignment of atoms along the attractive z-axis direction. A larger A value prevents 
such alignment which increases both kinetic as well as trap potential energy, hence increases the stability. In Fig. [ll], 
we display cm as a function of A, which separates the stable and unstable regions. Rigorous numerically calculations, 
on the other hand, find that condensate can still collapse even when A > A c . For example, we found collapse occurs 
when P > 5000 with A > 7, [c(£) > 2.0]. This difference maybe due to the choice of a simple Gaussian shaped 
variation function 



When the system is not in the Thomas-Fermi limit, we have to first solve for kq from equation 



4 A 2 = c(£ )[4 9 (ko) A 2 / M] + \ {K ° p ^F [1 " «o " c(£)(g(n ) ~ «g/(«o))] \ , 1 47 ! 



1/5 



and then using 

n 2 [1 - c(£)g( Ko )} - X 2 [l-c(£)/M] 



v M = P- 



X 2 



and vq = kqVzO to find the equilibrium widths. In this case, we find possibilities for one, two, three, and four or no 
roots of Ko depending on values of P, c(£), and A. The stability conditions for these roots, however, are similar to 
the TEA case as discussed before — there is at most only one stable solution. In Fig. [l^, we present kq, Vq, and v z q as 
functions of c{£). We see that vq decreases with increasing c{£), and the condensate collapses when vq goes to zero. 
Figure |l2|(c) displays the dependence of the condensate volume on £, where the volume is defined as the produce of its 
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effective widths in three separate dimensions. In terms of absolute values, numerically calculated volume (averaged 
width) differs from the variational result by a few times due to the multiplying effects of three widths. The shrinking 
volume increases the condensate density, which in turn can significantly increase the three body loss process, providing 
a potentially useful diagnosis tool M] . 

Figure [l^ shows cm as a function of P for different A values. We see that, for smaller P, a condensate with a small 
A can be more stable than condensates with larger A; while for larger P, a condensate with a larger A is always more 
stable than condensates with smaller A. This also confirmed by both variational (TEA) and numerical calculations. 

After we first submitted this paper, a paper on the same topic become available |2q| . We therefore compare our 
results with those from |26| in the next few figures. Figure [l4| displays the change of aspect ratio of the ground state 
for two extreme values of A = 0.1, 10 and for a small value of P = 5, and can be directly compared with the Fig. 3 
of [^6| . We note that essentially the same results were obtained as in |26j presumably because their neglect of the 
s-wave interaction simply corresponds to P = of our more general results. 

More interestingly, we show in Fig. [l5| for a larger value of P = 500. We find the aspect ratio now changes in the 
opposite direction with increasing dipole interaction strength. This reversal of aspect ratio with increasing values of 
P (due to increasing in atom numbers or s-wave scattering length a sc ) is due precisely to the physics of minimizing 
the total system free energy discussed earlier in the TFA. This phenomena was not observed in the simpler model 
of Ref. |26|. We also find that A c remains virtually independent of P at the same value as in the TFA: 5.170169, 
consistent with Ref. [E6|. 



B. Evolution of widths 

The evolution of condensate widths are found by numerically integrating Eq. (fllj). Assuming initially c{£) = 0, 
v(0) — Vq, v z (0) — v z0 , and v(Q) = v z (0) = 0, we can apply an electric field suddenly or slowly for t > 0. Under stable 
conditions for the widths, we choose to apply electric field suddenly. Otherwise the electric is increased linearly within 
a ramp-up time, and then kept constant. First, for the stable case, Figure |l^ shows condensate widths evolution up 
to t — 30 (it has been calculated up to t = 1000). When c(£ ) = 0, the condensate remains at its initially equilibrium 
state, and the widths are unchanged with time. When c(£) ^ 0, condensate widths oscillate with time, and prolonged 
numerical propagation indicates that we always have v > and v z > 0, i.e. the condensate is stable. 

We also see from these figures that the oscillation amplitudes increased with increasing c(£). Then finally at some 
stage, we could arrive at v < or v z < 0, signaling the condensate collapse. Figure |l7] indeed displays such cases 
when a linear ramp-on of the external electric field is applied. 



C. Small amplitude shape oscillations 



Once the equilibrium widths are found from numerically solving Eqs. (|44|), small amplitude oscillations can be 
studied by evaluating the matrix of the second order derivatives of the equivalent potential U(v x ,v y ,v z ) Eq. (|3£ 
We find that it takes the following symmetric form 



U n U u U 13 
U12 Uu U 13 

U13 U13 U33 



(48) 



where Uij — Uji due to nature of commuting derivative operations with different coordinates, and Un — U22 and 
U13 = U23 due to the cylindrical symmetry. We find the oscillation frequencies to be 



1 



^2,3= 



^/2 



U u + U12 + U33 ± x/t/fi + U( 2 + U^ 3 + 81^ + 2UuU 12 - 2UnU 3 3 - 2U 12 U 3 3 



1/2 



(49) 



where the expression for the matrix elements Uij are listed in Appendix Typical results and mode structure 
identifications J38| are given in Fig. [if We see that mode 1 and mode 3 are doubly degenerate when c(£ ) = 0. This 
is due to the additional symmetry Un — U33 and U12 = U13 for A = 1. 
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V. CONCLUSION 



In conclusion, we have performed a detailed study of trapped condensates with dipole interactions. We have devel- 
oped a general scheme for constructing effective pseudo-potentials for anisotropic interactions ]4lf| , which guarantees 
the agreement between the first order Born scattering amplitude from the pseudo-potential and the complete scatter- 
ing amplitude obtained from a multi-channel collision calculation. Our theory has been applied to the study of induced 
electric dipole interactions and can also be directly extended to magnetic dipole interactions as well as permanent 
electric dipole interactions of trapped molecules p7| . 

Finally we provide a reality check for prospects of experimental observations of the electric field induced interaction 
effects. Though the required fields arc relatively high, there are evidences they can be created with current laboratory 
technology. In Ref. [Q fields of upto 2 x 10 5 (V/cm) were used to slow a molecular beam. Gould |^3| used fields upto 
4.6 x 10 5 (V/cm) in the measurement of atomic tensor polarizability, while Marrus et al. [Q reported fields upto 10 6 
(V/cm). What is perhaps most encouraging is a recent experiment for cooling molecule beams with time-dependent 
(adiabatic from the view point of atomic internal dynamics) fields of upto 10 7 (V/cm) [ p5[ . We also note that at the 
high fields being discussed in this paper, the tunneling ionization of atoms remain infinitesimally small |4q| . 

This work is supported by the U.S. Office of Naval Research grant No. 14-97-1-0633 and by the NSF grant No. 
PHY-9722410. The computation of this work is partially supported by NSF through a grant for the ITAMP at 
Harvard University and Smithsonian Astrophysical Observatory. 



APPENDIX A: U-MATRIX ELEMENTS 



After tedious calculations, we find that 



P 



VqV z0 



5ttc(£) 



24(«2o-«g)3 



54v 2 v 4 z0 + 16v% - 9(llv 2 + 4v 2 z0 )v 4 H(v /v z0 )) 



(Al) 



U: 



= A2 + -r 



P 



5ttc(£) 



3(^ 2 o-«o 2 ) 3 



(4« 6 - 12vy z0 + 51vy z0 + 2v% - 9(v 2 + Av 2 z0 )vy z0 H(v /v z0 )) 



(A2) 



u- 



p 

12 — —a 

VqV z0 



5ttc(£) 



l&vt + 51 V y z0 - 30vy z0 + 8v% - 45v 6 H(v /v z0 )) 



(A3) 



and 



U 



P 



13 



,,3„,2 



5ttc(£) 



6(?4 - v o) 3 



(4« 6 - 36vy z0 - 15v 2 y zQ + 2v% + 45vy z0 H(v /v z0 )) 



(A4) 
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FIG. 2. Computation of kernel 
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FIG. 3. Typical behavior of K.(p' , p, 0) (a) and the corresponding integrated kernel over the wave function grid (b). The solid 
(dashed) lines are for p — 1.0 (3.0). 
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FIG. 4. cm at the highest field value of condensate collapsing for A = 1. Numerically computed values are in circles connected 
by a dotted line. The solid line is from the variational calculation. 
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FIG. 6. Ground state wave functions with (solid line) and without (dash-dotted line) TFA. P — 5000, A = 1 and 
c(£ ) = 0.7, 0.6, 0.2 in descending order of wave function values at the origin. 
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FIG. 8. Condensates (pancake shaped in the left and cigar shaped to the right) always expands along the direction of an 
externally applied electric field and shrinks along the repulsive radial directions. Darker ellipses are for higher electric fields. 
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FIG. 11. cm as a function of A with TFA. 
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FIG. 14. Field dependence of ko and condensate widths (in insets) for P — 5, A = 0.1 (a) and A = 10 (b). When A < A c . 
volume of condensate decreases with c(£). It increases with c(£) when A > A c , causing condensate to be always stable. 
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FIG. 16. Evolution of condensate widths for P = 10, A = 1, v(Q) = w z (0) = 1.63359 at c(S) = 0.0 (solid line), c{£) 
(dashed line), and c(£) — 0.8 (dash-dotted line). 
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FIG. 17. Evolution of condensate widths at c{£ ) = 1.0 (> cju = 0.9989), for electric field ramp-up time T = 5 (solid line), 
10 (dashed line), and 20 (dash-dotted lines). Other parameters are P — 10, A = 1, and v(0) = v z (0) = 1.63359. 




FIG. 18. Electric field dependence of the shape oscillation frequencies. Other parameters are P = 1 and A = 1. 
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